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ABSTRACT 


This thesis addresses three deficiencies in the NPSNET simulated world. First, 
although a full set of algorithms have been defined for Dead Reckoning (DR) entities in a 
distributed simulation, NPSNET only implements a few simple linear algorithms. Second, 
NPSNET lacks a set of physically-based models for munition trajectories having, currently, 
only algorithms for the bullet and bomb. Third, NPSNET lacks physically-based models 
for engine power curves using, instead, a simple linear approximation. 

The purpose of this thesis work is to implement an object-oriented programming 
toolkit which corrects these deficiencies. The code, in C++, utilizes class hierarchies. The 
toolkit implements the nine class hierarchies of DR algorithms described by the Advanced 
Research Projects Agency for the Distributed Interactive Simulation standard. The toolkit 
also provides treatment of a physically-based class hierarchy for munitions trajectories. In 
addition, a physically-based, mathematical model for the engines class was implemented. 

As a result, a set of DR algorithms have been built to predict the position of 
simulated entities in all cases. The munitions class implements trajectories for a variety of 
projectiles. With this arsenal, future versions of NPSNET will be more realistic. The engine 
class, with new mathematical models, far more realistically represents engine behaviors 
than the current linear approximation. In summation, the implementation of this toolkit 


dovetails very well with the needs of NPSNET, and will be integrated into future releases. 
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I. INTRODUCTION 


A. BACKGROUND 


Conventional training and preparation for future military operations is an expensive 
and time-consuming process involving many persons and resources. With the current 
budgetary constraints and downsizing of the armed forces, ways of doing business must be 
revised, but in a way that does not compromise preparedness. Fortunately, with advances 
in Computer Science, military operations from battle to replenishment can now be 
simulated on computers, and personnel can complete many phases of training in simulated 
battlefields, thus mitigating the expense of engaging in the real-world environment. 

One such simulated battlefield, NPSNET [PRAT93], has been developed by the 
Computer Science Department of the Naval Postgraduate School (NPS). NPSNET is a low- 
cost, distributed, interactive, virtual 3D world for simulation and training. Personnel in the 
simulated environment interact with each other as if they were on a battlefield. They can 
virtually fly an airplane, ride a tank, walk, run, and shoot at each other. 

NPSNET is an evolving system, a test-bed for new research as well as a growing 
amalgam of proven ideas. Recently, numerous projects have improved and enhanced 
NPSNET. Such works encompass: Configuration Management, Real-time Scene 
Management, Terrain Development, Insertion of the Human, Human-Computer Interface, 
and Networking. However, NPSNET lacks true physical representation of moving bodies. 
Only a few simple cases have been implemented such as linear Dead Reckoning (DR), and 
simplified munitions trajectories for the bullet and bomb. Improving the physically-based 


aspects of NPSNET is the purpose of this thesis work. 


B. PROBLEM STATEMENT 


This thesis addresses and resolves three specific deficiencies in the physically- 
based representations of NPSNET. The first is dead reckoning (DR). NPSNET uses the 


Distributed Interactive Simulation (DIS) protocol, a standard for networked simulators, to 








Communicate between entities--aircraft, vehicles, munitions, etc. -- in its environment 
[IST93]. Every workstation participating in a simulation, each running NPSNET, shares 
the data for the entities it controls, thus allowing each NPSNET to display a consistent view 
of the simulated world to its user. Each entity in the environment sends to the network its 
velocity, position, and orientation no less frequently than every five seconds. Other objects 
receive and render that information within their own workstation. For efficiency and to 
conserve network bandwidth, data for an entity is not sent out continually. To keep an 
updated picture of networked entities, each workstation dead reckons (predicts) remote 
entities’ velocity, position, and orientation based on its last update. Specific dead reckoning 
algorithms constitute part of the DIS standard. Currently, NPSNET implements only a few 
simple linear world coordinate cases. A full implementation of DR for all vehicle cases in 
both world and body coordinates is required. 

The second deficiency regards simulation of munitions. In a simulated battle, 
continuous simulation of munitions is absolutely essential. Currently in NPSNET, only a 
few simple cases of a bullet and a bomb are implemented. A more definitive physically- 
based implementation of munitions is required to more closely reflect real-world munitions 
behavior. 

The third deficiency is engine simulation. To make the virtual world realistic and 
plausible to the user, moving objects must reflect the behavior of their real-world 
counterparts as much as possible. Currently, NPSNET lacks any engine type. Therefore, a 


basic implementation of engines is needed, too. 


C. APPROACH 


For DR implementation, we use the study of the Advanced Research Projects 
Agency (ARPA) for DR algorithms published as part of DIS 2.0.3 [TOWE94]. In this 
study, the equations of motion for the DR algorithms are divided into nine groups. The first 
group is for a static case. Another four groups are for world coordinates, and the remainder 
are for body coordinates. Using the C++ object-oriented programming methodology, this 


thesis work implements all nine groups into a class hierarchy. 














The munition implementation is based on Newtonian Mechanics. Most moving 
objects in the real world are governed by gravitational forces. Therefore, kinematic 
equations with constant downward acceleration are utilized for bullet, bomb, ballistic, and 
artillery objects. A physically-based model of a system with lost mass is used for rockets 
and missiles. Air drag 1s omitted. 

For engine implementation, we built an abstract input and output system with 
general quadratic, logarithmic, and exponential models of acceleration. (A full physically- 
based model is beyond the scope of this thesis.) 

All of the implemented code is written in C++ with class hierarchies for portability 
and reusability. With its class structure and inheritance, C++ minimizes code development 


as well as making the program easier to understand. 


D. ORGANIZATION 


This first chapter introduced NPSNET and discussed the problems addressed by 
this thesis. It also described the approach taken to solve each problem.The second chapter 
covers the background and previous works impinging on the problems to be solved. The 
third chapter describes in detail the implementation of the DR standard. The fourth chapter 
recounts the formulas and equations which govern the munition classes. The fifth chapter 
describes the implementation of the engine mathematical models and the derived formulas. 


The sixth chapter contains the conclusion and presents necessary future work. 











II. BACKGROUND AND PREVIOUS WORK 


A. PREVIOUS WORKS 


1. Dead Reckoning (DR) 
Currently, NPSNET implements only a few simple cases of DR where there 1s no 
rotation involved. There are the static case, the two cases of world-coordinates, and the 


default case for the rest of the situations, see Figure 1. 


moveDR 


Case 3 


Case 1 r=0, v=x, a=x 
r=0, v=0, a=0 Sener se Ue 


Case 2 | Default 
Case 


r=0, v=x, a=0 
2 , vV=x 


r = rotation a = acceleration 
v = velocity x = Some number 





Figure 1: Current NPSNET DR IMPLEMENTATION 


The first algorithm applies to a static object when its rotation, velocity, and 
acceleration are all zero. The second algorithm applies when object rotation and 
acceleration are zero but velocity is not. The third algorithm applies when object rotation 
is zero and velocity and acceleration are not. The default case applies when object velocity 


is a constant, so that the object moves in a straight line. These four DR cases were 








implemented in C++ as special cases in the subclass moveDR of the Vehicle Class in 
NPSNET. Kinematic equations with constant acceleration were used to implement those 


Cases: 


v, =v, +a;t Eq 1 


E Qo 
x, = x, +v,t+0.5a,t? Eq 2 
Where: 


v, = respective velocity in each direction 


v, = respective initial velocity in each direction 


oO 


x, = respective position in each direction 
x, = respective initial position in each direction 
a, = respective acceleration of each direction 


t = time 


One of the purposes of DIS, a set of standards developed by ARPA and industry, is 
to provide a specification to be used by government agencies and engineers to build 
interoperable networked simulation systems [IST93]. DIS defined a need for DR, and 
ARPA released its DR algorithms study on February 7, 1994 [TOWE94]. In this study, with 
C++ methodology in mind, the algorithms are divided into a general static case and two 
large groups, one for the world-coordinates, and the other for body-coordinates. Inside each 
large group are four subgroups, each addressing different cases of the moving object in its 


respective coordinates, see Figure 2. 








World coord. Body coord. 
Figure 2: DR algorithms of ARPA’s study 





2. Munition 
In the current version of NPSNET, only a few simple cases of munitions are 
implemented, such as a bullet, a bomb, and a simple non-loss-mass missile, see Figure 3. 


Equations similar to Eq1 and Eq2 above are used to implement those munitions. 


Simple 
Missile 


Figure 3: Current NPSNET munitions implementation 





Newtonian Mechanics governs the behavior of objects flying in the air in a constant 
gravitational field. Equations that describe motion and kinematics have been developed 
ever since Newton presented his greatest discovery (F = ma). Many first-year math and 


physics books describe these equations. [MART84][HALL78] 





Another relevant work is the NATO study of ballistic missiles [NATO]. It is the 
definitive study of the trajectory of a flying mass in the air. It describes all factors 
influencing the behavior of flying objects, encompassing drag force, lift force, Magnus 
force, and Coriolis force. Most of these equations exceed the scope of this thesis and are 


unnecessary for NPSNET, though it could apply to later enhancement work. 


3. Engines 

Engine models are not currently incorporated into NPSNET. Presently, NPSNET 
uses a simple linear approximation between the user input to a throttle and the speed of the 
entity. However, the real world of vehicles encompasses a variety of engine models, such 
as motors, gas turbine, piston, and rocket engines. All currently existing engine models 
exceed the scope of this thesis, except for the rocket engine. With a few revisions, it follows 


the rocket munition model and so is implemented in this thesis work. 
B. SUMMARY 


In summation, NPSNET lacks complete physically-based models for DR, 
munitions, and engines, for accurately representing objects in its virtual world. All DR and 
munition models are readily available, and can be implemented. Except for the rocket 
engine which can be implemented directly from rocket models, the remaining engine 
models need to be abstracted as an input and output system and tested until it seems 


plausible in the virtual world. 








HI. DEAD RECKONING 


A. OVERVIEW 


What is dead reckoning? According to the Random House dictionary, dead 
reckoning is the “calculation of one’s position on the basis of distance run on various 
headings since the last precisely observed position.” Our usage of the word dead reckoning 
in this chapter is exactly the same as the definition. 

To allow multiple users in a simulation, we run NPSNET on different workstations 
which are all connected by the same network. Each workstation presents a consistent view 
of the entire simulated world, rendering both the local objects it controls and the networked 
objects controlled by the other simulators. Each simulator transmits the position, velocity, 
and orientation for its own local objects, and reads the network for information on remote 
objects. This state change data is sent on the network no less frequently than every five 
seconds in a packet called, in DIS terminology, a Program Data Unit (PDU). But five 
seconds is a long time when the user is watching an entity which is animated on the screen. 
The entity cannot jump to a new location every so many seconds and satisfy the appearance 
of smooth motion. In-between times, each workstation needs to dead reckon each remote 
object’s position, velocity, and orientation and move it smoothly along its probable course, 


Figure 4. 





Figure 4: PDU and DR time 





In endeavoring to provide a complete set of DR functionality, this thesis utilizes the 
algorithms developed by ARPA for both world and body coordinates. These are 


implemented using appropriate C++ classes. 


B. DESIGN 


STATIC 
Dead Reckoning 


Algorithm I 


RMI) — CDRME> 


WORLD COORDINATES BODY COORDINATES 





Figure 5: DR classes 


In object-oriented languages like C++, inheritance is when a child class inherits the 
variables and functions of a parent class. Thus the child class is automatically like its parent 
except for the specific ways in which it needs to be different, which makes design and 
coding clean and efficient. In Figure 5, each of the nine ellipses represents a class of DR 
algorithm. The parent class is the static class. The four classes on the left side of the 
structure belong to the world coordinates class, the four on the right belong to the body 
coordinates. 

The structure in Figure 5 varies a bit from ARPA’s DR study, see Figure 2. As 
represented in Figure 2, dead reckoning method four (DRM4) can be derived from either 
DRM3 or DRMS; likewise for DRM8, which can be derived from DRM7 or DRM9. But, 
in Figure 5, DRM4 is only derived from DRM3, and DRM8 is only derived from DRM7. 
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The reason for the variance is to reduce the complications arising from muitiple inherited 
variables and functions. It is far more straightforward to keep track of ali the variables and 
functions if each child class derives from a single parent class. 

Table 1 lists each DR algorithm in the set. A key to the abbreviations: 


B = body coordinate 
F = fixed 

P = position 

R = rotation 

V = velocity 

W = word coordinate 


For example, the fourth row indicates that algorithm DRM4 has rotation, velocity, 
and belongs to the world coordinates subclass. Its roli, pitch, yaw (orientation), velocity, 


and acceleration are equal to some number. 
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Table 1. DR Algorithms 
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Furthermore, in the representation DRM(X,Y,Z), X represents for the rotation of 


the entity, Y represent for position, and Z represent for coordinate. 
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Cc. IMPLEMENTATION 


1. Static class (DRM1) 

This is the parent class for all other DR classes. Its rotation, velocity, and 
acceleration are all zero, so its position is fixed. To use it, we simply instantiate an object 
and set its position to given values. In this case, we always get back the same initial value, 


when we call for updated position. 
2. World-coordinates (DRM2-5) classes 


The input to all the algorithms are: 


X ol t | = position vector in world coordinates at initial time. 


velocity vector in world coordinates at initial time. 


< 

= 

ae 
II 


A,, = acceleration vector in world coordinates. 
(Cy (t,) 8 (¢,) 0 (¢,)) = Euler angles at initial time. 
@ = (© ,@ o,)F = (01.0, 7 = angular velocity in body coordinates. 


At = time increment for dead reckoning step. 


The outputs are: 
Xl % + At) = position vector in world coordinates after time delta. 


Vee, Pas At | = velocity vector in world coordinates after time delta. 


(v( t + Ar},6( as Ar),o t + At) | = Euler angles after time delta. 


(1) DRM2 
In this algorithm, r = 0, v = x, and a=0. We use the Eq3 and Eq4 to compute 


the object velocity and position at time ¢,+Ar. Euler angles at time delta 
(v(t, + Ar},0( 4, +A0).0 1, +A) are the same as its’ initial time 


( v( t » 9 t, 9 5) because there is no rotation involved. 
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V t + Ar) = V t )+ArA Eq3 
Wi O Wi O WwW 


x, [ t + Ar) a X(t) + atv,,{ ty) + SAA Eq 4 
(2) DRM3 
In this algorithm, r = x, v = x, anda =0. We can get velocity and position at 
time t,+At by Eq3 and Eq4, and the steps to get Euler angles at time ¢, + At are as 
follows: 
Step 1: Use the initial Euler angles (y (t,) ,O8(t,) ,0 (¢,) ) to compute the rotations 


matrix U (b(t,) |W) which takes world coordinates into body coordinates at time t, . The 


formula for U (b (t,) |“) 1S 


e(w) (8) s(w)c(8) -s (8) 
U(4(1,} pw) = |-scw)c(@) +5(6)5(8)eCw) e(O)e(W) +5(6)s(@)s(W)_s(@)e(®)| FAS 
s()s(W) +¢(9)s(@)c(y) —s() eC) +€ (0) 6(8) s(W) ¢ (0) ¢ (8) 
Where c = cosine and s = sine. 


Step 2: Compute the rotation matrix U (b(t, + At) [5 (t,)) which takes body 


coordinates at time ¢, into body coordinates at time ¢, + At. 


U(o(1,+41)10( #,)) = (= cos I@lAD) oT 4 COs (lola 1- SP GEID Eg 6 


jeo|? 
Where: 


The matrix I is the 3x3 identity matrix having 1’s on the main diagonal and 0’s off 


of the main diagonal. The matrix ww! is 
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©, 0, 0,0 W, 


2 3 
— 
30, 0, 0303 
and the magnitude of @ is 
lo] = (0? + 2 + 03) 1/7 Eq9 


Step 3: Compute the rotation matrix U (b(t, + At) \“) which takes world 


coordinates into body coordinates at time t, + At. 


U(b(t,+At)|w) = U(b(t, + At) |b(t,)) -U (b (t,) [w) Eq 10 


Step 4: Compute the Euler angles (wrt, + At) ,0(t, + At) 0 (1, + At)) after the 


time increment using the rotation matrix U(b(t,+ At) \w) . The rotation matrix 


U (b(t, + At) |”) has the form: 


O11 Uy. U 43 
U (b(t, + At) |w) = |U,, Uy Us, Eq 11 
Us, Us U 33 
Where the entries U.. are functions of the Euler angles 


Ly 


(y(t, + At) 6 (t, + At), (t, + At) ) as shown in the matrix displayed in Step 3 above. 


Step 5: The Euler angle @ is recovered via 
—T 1 
Q = —asin (U,,) € | Eq 12 
If Us and Ui. do not both equal zero, the Euler angle y is recovered via 


w = Arg (U,,,U,,) € [-1,7] Eq 13 
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Where the function Arg(x,y) is an argument of the complex number x + yi, the 
particular branch of the argument chosen to lie in the interval [—1,7] . 


If z7.. and zy,. do not both vanish, the Euler angle 9 1s recovered via 
23 33 
o = Arg (U,,,U43) Eq 14 


The exceptional cases occur when 9 = 1/2 or —1/2. In that case, the Euler angles 
w and o are recovered via 
y = 0 Eq 15 


o = Arg (U,,,-U3,) Eq 16 


(1) DRM4 
In this algorithm, r = x, v = x, and a = x. Position, velocity, and orientation 


are computed much the same as DRM3 but with the addition of acceleration A,, . 


(3) DRMS5 
In this algorithm, r = 0, v = x, and a= x. We can use Eq3 and Eq4 above to 
compute its position and velocity. DRMS lacks rotation so we do not require equations Eq5 


through Eq16. Its orientation is fixed. 


3. Body-coordinates (DRM6-9) classes 
The input to all the algorithms are: 

X,, (tQ) 

Vy (t,) 


A, = time derivative of V, 


position vector in world coordinates at initial time. 


velocity vector in body coordinates at initial time. 


(y (t,) 8 (¢,) .6(¢,)) = Euler angles at initial time. 
QO = (o,.0,00, \P = (005.03 |! = angular velocity in body coordinates. 


At = time increment for dead reckoning step. 
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The outputs are: 


xX wl Ls Ar = position vector in world coordinates after time delta. 


Al 


( ¥{ al Ar),6{ ¢, + Ar),o( t+ At) = Euler angles after time delta. 


t+ At | = velocity vector in body coordinates after time delta. 


(4) DRM6 


In this algorithm, r = 0, v = x, and a = 0. We can compute the velocity and 


position at ¢, + At with Eq17 and 18, ignoring second terms on the right of both equations 
because they are zero. The Euler angles at time t, + At is the same as its initial condition 


(yw (t,) 8 (¢,) 0 (¢,)) , because there is no rotation involved. 
Vi 40 + Ar} = Vat | + At, Eq 17 
X[ 5 + At) = X,( 45) + u(mie(«,)\ 9 av, ¢,) + [ofce™@dca, | Eq 18 


foterta, = WlAr—sin (lalAd) gr, sin (lala ; , 1— cos (QJAD 9 54 19 
0 |co|> lo lal? 


slol2ar — cos (|@|Az) — |e Atsin (lalAz) + 1 


T 
lol? M@* + 


[jy tetra = 


cos (|@|Ar) + os (|@| At) — 1, , sin (|co| Ar) cea (|ao] At) Q £420 


Computing |@| and Q is the same as in Eq9 and Eq7, and U (12( , }] is the 


transposition of the matrix U (b(t,) |W) above. See Eq)5. 
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The following steps are used to calculate Euler angles 


(y(t, + At) 0 (¢,+ At) .o(t,+At)) attime t+ At. 

Step 1: Compute the rotation matrix U (b(t, + Af) | (t,)) which takes body 
coordinates at time f¢, into body coordinates at time t, + At, the same as Eq 6. 

Step 2: Compute the rotation matrix U (b(t, + At) |) which takes world 
coordinates into body coordinates at time t, + At as Eq 10. 

Step 3: Compute the Euler angles (y(t, + At) ,6(t, + At) ,0(¢,+Ad)) after the 
time increment using the rotation matrix U (b(t, + At) \w) , same as the Step 4 in the 


world coordinates. 


(5) DRM7 
In this algorithm, r = x, v = x, a= 0. We can compute the velocity, position, 


and Euler angles at t, + At with Eq3, Eq4, and Steps | to 3. 


(1) DRM8& 
The steps to compute the outputs for this algorithm, where r = x, v = x, and 


a =x, are the same as the steps of DRM7 with the exception that a is not zero. 


(6) DRM9 
In this algorithm, r = 0, v =x, and a =x. The steps to get the outputs are the 
same as DRM6 with a not equal to that zero. The only thing changes in this algorithm is 


that the second term on the right of Eq17 is not zero. 
D. RESULTS 


With the organization afforded by the C++ class hierarchy, the coding and testing 
procedures for DR classes are segregated into separate, discrete units. The debugging 
process for each implemented algorithm was very simple. For each case, the code was 


tested, as suggested in Figure 4, as follows (see code in Appendix A for details): 
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a. Instantiate an object of that case. | 

b. Call procedure update_DR to set the object to its “real” values (as the information 
is obtained from incoming PDUs). 

c. Intermittently call move DR to compute the predicted position, velocity, and 
angular position of the object. 

d. Repeat step c. at the desired interval until the next PDU arrives with real data for 
the object, then repeat step b. 

The code is extremely reliable, and its performance very quick, even when tested 
on a low-end Intel 80486-based PC. To illustrate test results, the following Tables 2 and 3 
list the input and output tested cases for algorithms DRM4 and DRM8. They are the two 


most complex cases in the DR classes. 


[Initial Conditions [x | 




















x y 

ion miseoreeo] 80 | 05 _| 90 
Rotation Rate -degiee| 0 | 6 | 70 
- FinatResus | | | 


Table 2. Test of DR algorithm 4, Delta time is 0.5 sec 


0 
Acceleration - m/sec*sec 
5 
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Table 3. Test of DR algorithm 8, delta time is 0.5 sec 


Following Table 4 is a time performance data comparison between this thesis’s 
code implementation in C++ and ARPA’s study code implementation in C. Those 
implementations were run on different platforms to test the consistency and time usage of 


the code. 


SGI Indigo 2 | SGI Indigo 2 
1S0MHz 1S50MHz 
Optimized 





Table 4. Performance Data 
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Observe that the trend of the times spent on computation were similar in all cases 
and different platforms. Without optimization in compiling, this thesis implementation ran 
fast in the simple cases of algorithms 1, 2, 5, 6, and 9, but slow in the complicated cases of 
3, 4, 7, and 8. Using compilation optimization (with the flag “-O3 -mips2”), and 
considering the lower overhead resulting from putting each class into its own file, most of 
the cases ran faster than the ARPA’s study code. However, in the complex case of 
algorithm 8, our code ran a little slower due to redundant matrix multiplication in 
computing angular position. The speed of the 486DLC 33Mhz CPU is about three times 
slower than the Intel 486DX2 66MHz CPU. 

In summation, the code for each of the DR algorithms works as designed, and has 
been implemented in the current NPSNET. NPSNET thus takes an important step forward 
in DIS-compliance. NPSNET can move networked entities from disparate simulators 


following any of the DR algorithms those entities specify. 
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IV. MUNITIONS 


A. OVERVIEW 


One of the goals of simulation is to make the virtual world resemble the real world 
as Closely as possible. This chapter describes an extensive munitions class, built so that 
NPSNET programmers can enhance and improve the physically-based representation of 
munitions in NPSNET. The class incorporates the mathematics that will allow various 


munitions to follow realistic trajectories in the simulated world. 
The kinematic equations v = v,+at and x = x,+v,t+0.5at* calculate velocity 


and position for most of the munitions, thus forming the core of the base class. Orientation 
of the munition is derived from its current velocity vector. Implementation of the 
subclasses of the base class is straightforward from those equations, except for the rocket 


model which involves a second order differential equation. 
B. DESIGN 


The munition class is the base class containing all the common variables and 
functions for its subclasses, see Figure 6. The equations that govern the velocity and 
position of the bullet, bomb, and ballistic subclasses are much the same. In the artillery 
class, though, we find the artillery’s necessary initial velocity first, then it uses the same 
equation as its parent class--bomb/ballistic--to find its velocity and position at a later time. 
The rocket subclass is totally different from its peers. Its velocity and position are governed 
by a different set of complex equations. The approach to compute the velocity and position 
of the ballistic missile class is the same as for the artillery class. That is, its initial conditions 
are determined before using the rocket’s velocity and position functions to find its velocity 
and position at a later time. 

The guided missile seeks and destroys its target. It is a subclass of the base class. 
Incorporated with the engine class, its trajectory can be further refined by factoring in the 


power curve as it is affected by user input. As examples, in one case the missile may couple 
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Bomb/ballistic Guided Missile > Missile 


Ballistic missile 





Figure 6: Munition class structure 


with a constant speed curve to chase after the target, while in another it couples with the 


quadratic engine speed equation. 


C. IMPLEMENTATION 


1. Orientation 


The orientation of the munition is converted from its velocity vector, a function 
residing in the munition base class. The function takes a velocity vector and returns angular 
position in degrees, the values commonly used by graphics software. This function checks 
all eight cases of given (x, y, z) values, and then used inverse tangent or direct assignment 
to assign appropriate angular position for each direction. 


For example, if the input vector is (x, y, z) then: 


6 =rtod- atan~ Eq 21 
x y 

8 =r1tod- atan~ Eq 22 
y x 
Z 

@ =r to_d- atan- Eq 23 
Z x 


Where r_to_d is the factor to change from radian to degrees. If the input vector is 


(O, y, z) then: 
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x 
6 =rtod- atan~ Eq 25 
y Z 
Z 
6 =r to_d- atan- Eq 26 
= y 


2. Bullet 


A bullet is a small projectile fired from small arms. To simulate its path, we use 
equations Eq 27 through Eq 32. The gravity constant is set to one or zero to implement the 


bullet with or without the effect of gravity, see Figure 7 and equations Eq 27 through Eq 32. 


bullet path without gravity 


™~ 
bullet path with gravity 





Figure 7: Bullet path 


Inputs: 

oF initial velocity in the respective direction. 
X, initial position in the respective direction. 
Outputs: 


v,(t) velocity at time t in the respective direction. 
x,(t) position at time t in the respective direction. 


@;(t) angular position at time t in the respective direction. 


The following equations Eq 27 through Eq 29 describe the velocity of the bullet. 


Velocity in the x and y axes of the horizontal plane will not change from the initial state 
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because there is no significant force acting on it. We omit the air drag and wind effects, 


which are negligible. Gravitational force acts upon the z direction. 


x x 
O 
voe=vy Eg 28 
J Le 
Vie —gt- gravity Eq 29 
O 


The following equations Eq 30 through Eq 32 describe the position of the bullet. 
Position, at any time, is based on initial position and initial velocity. In the z direction, the 


bullet drops due to gravitational force. 


Xx =x. +v_t Eg 30 
x x x 
O O 
x =x +VvV it Eq 31 
y yy 2h 
5 ite, mn t—0.5e@t* - gravity Eq 32 
O O 


3. Bomb 

Typically, a bomb is dropped from a flying object. Due to gravity, it descends from 
its initial position and accelerates toward the surface of the earth. Its path, normally, is a 
parabola, Figure 8. The equations governing its velocity and position are much the same as 
the equations for the bullet, Eq 27 through Eq 32. Inputs and outputs are the same as for the 


bullet class above. 


Together with Eq 21 and Eq 22, the following equation describes the velocity of the 
bomb. 











initial velocity 


initial position 


final position 





Figure 8: Path of a bomb 


Y= Ue 0b Eq 33 


Together with Eq 30and Eq 31, the following equation describes the position of the 
bomb. 


= 4 Z 
x= x, + vs O.5 gt Eq 34 


These equations are the same as Eq 27 through Eq 32. The difference is the 
“gravity” constant of Eq 29 and Eq 32 is omitted, since bombs are always dropped from 


high above the surface where gravitational force always applies. 


4. Ballistic 


Ballistic behavior applies to projectiles, like an artillery round, which are shot not 
directly at the target, but along a parabola whose path is calculated to culminate at the target 


location, Figure 9. Inputs and outputs are the same as for the bullet and bomb classes above. 
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final position 


initial position 





Figure 9: Path of the ballistic 
Equations that govern the behavior of ballistic projectiles are much the same as 
equations for the bomb. The only difference between the two classes are the initial 
conditions, that is, the ballistic projectile normally originates from the ground and the bomb 


from the air. 


a. Artillery 

The path of the artillery round is the same as for the general ballistic round. 
The difference is in the given conditions. For the ballistic, we are given the initial velocity, 
and for the artillery we are given the artillery and the target positions. 

Inputs: 

P1(x1, yl, z1) initial position of the artillery. 

P2(x2, y2, Z2) initial position of the target. 


Outputs: the same as for the bullet. 


In the following case, we only know two positions: P1 and P2. We do not 
know the artillery round’s initial velocity, see Figure 8. If we did know the artillery’s initial 
velocity, the way to find the position at time t, P(t), would be the same as for the bomb or 
ballistic. So the approach is to find the initial velocity first. Once we have the initial 
condition, the equations that govern the behavior of the velocity and the position at a later 


time are the same as the equations for the bomb or ballistic. 
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Figure 10: Artillery’s round path 


The following equations Eq 35 through Eq 37 are used to find the required 


nozzle velocity of the artillery round. 








X—X, 
y = Eq 35 
XQ t 

ye 
y = Eq 36 
Yo t : 


ea ie 0.5¢t7 
uy = a as Eq 37 
In Eq 35 through Eq 37, the unknown variable on the right is the time, t, so 
the programmer can code the implementation to query the user for that value. In Eq 37, the 
top right part takes the absolute value to ensure that the velocity in the z direction is positive 


(going up). 
5. Rocket 


A rocket is a type of munition that burns fuel, shooting it rearwards to gain velocity 
and altitude, hence distance (Figure 11). At a certain burn-out point, the rocket will have 
expended its fuel. From that point, it enters free-fall as a ballistic projectile. Its subsequent 


path is based on the velocity and position gained in the burn phase. Therefore, there are two 
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free falkregion 





Figure 11: Rocket’s path 
set of equations governing the rocket’s velocity and position, the first applying to the burn 


phase, the second to the free-fall phase. 


Inputs: all the variables on the right side of the following velocity and position 
equations. 


Outputs: the same as before: velocity, position, and angular position. 


Equation 38 states that the rate of change of the momentum (P = mv) is equal to the 


sum of the external forces that act on the object. 


ad 
P 7 Poet Eq 38 


Based on Eq 38, we can derive the equation for the motion of an object with lost 
mass, Eq 39, where m is mass, V is velocity, u is the velocity of the fuel, t is time, and F is 
the external force. Eq 39 states that mass multiplied by the rate of change of velocity, minus 
velocity multiplied by the rate of change of mass, equals the external force that applies to 


the object. 
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d dm ) 
mV —b =F, Eq 39 


From Eq 39, and if we assume that a r where r, rate of lost mass, 1s a constant, 


we can derive the velocity and position equations. The derivation can be found in many 


math and physics books [MART84][HALL78][BOYC77]. 


Eq 40 through Eq 42 following are the velocity equations for the rocket. 


v—-v =U nf 1 | Eq 40 
vA m 
O O 
rt 
v.—-vV. =U nf 1-5] Eq 41 
x x x m 
O O 
rt 
vy —v. =u In| 1-— Eg 42 
ye y mM, 


Eq 43 through Eq 45 following are the equations for the position of the rocket. 











ae rt rt ] 
ek <9. {-[-ame Jt tig | a8 Eq 43 
zg Z r m m 2 
O O O O 
Yo rt rt 
xX -X. =v t— [-[ me [ae |] Eq 44 
x “Xx x r m 
O O O O 
6 rt rt 
x-x =v t-2 {1 ] —-— hes ]-—— Eq 45 
y yy i y b r Le m 
Where: 
= respective velocity in each direction 
R. = respective position in each direction 
i. = respective fuel velocity in each direction 
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m = total mass of the rocket 


r = rate of lost mass 
= gravitational constant, 9.81m/s2 


§ 

F zt = respective initial velocity and position 
0 0 

t = time 


If we assume that fuel mass will burn-out at a certain time, as in Eq 46, we can 


derive the burn-out velocity and the burn-out position as in equations Eq 47 through Eq 52. 


m =m —m = rt Eq 46 
fuel O rocket 
Where: 
m uel = fuel mass 
m = mass of the bare rocket 
rocket 


The following equations Eq 47 through Eq 49 give the burn-out velocity of the 


rocket. 
m m 
v—-v =-uU of 1+ |e | Eq 47 
a *Z Z m r 
O r 
ei 
vy —-v. =-u Inf 1+— Eq 48 
x x m 
O r 
a ee “f 
vy —v. =-—u In| 1+ Eq 49 
yw mM 
The following equations Eq 50 through Eq 52 give the burn-out position of the 
rocket. 
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Pas, sa Pee of -(Ze][s—a(-Z2)H} 4) Eq 50 
is Z r m m pean 
O O O r 
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ae : l-| — fin 4 Eq 51 
- i r mM, m 
um Mm, mM 
x ES ea | ee ee Eq 52 
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a. Ballistic missile 

The inter-continental ballistic missile (ICBM) is a type of weapon which 
uses rocket thrust and has the capability of flying directly to a given target position, Figure 
12. The equations governing its velocity and position are the same as for the artillery class. 
We determine the ICBM initial conditions. Then we use the parent class velocity and 


position equations to find the velocity and position at a later time. 


burn-oxt point 


free falkregion 


Given: pl and p2 





Figure 12: Rocket shot at predetermined target. 
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6. Guided missile 





Figure 13: Missile chasing target 


A guided missile is a projectile with the capability of chasing and homing 
in on a moving target, Figure 13. Our general solution for the moving target will apply to 
the fixed target case, too. In the moving target case, the missile must change its velocity 
vector every delta time, At = t2 - tl. At the end of each delta time, we know the position 
and velocity of both the missile and target. We use that as the initial conditions for the next 
time interval. We repeat the process until the missile is within a detonation range from the 
target. 


The following equation Eq 53 computes the new missile velocity for every At. 


Preset wile 
VD. ap Speed Eq 53 
missile |? = ee | 
target “ missile 
Where: 
ie arget = ‘atget position 
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> . ., = new missile velocity 
missile 

wieie russe position 

Speed = speed of the missile from previous delta time 


Following is the position equation for the guided missile 
et + y; At + 0.54, (At) 2 Eq 54 
Where i = is the respective direction, x, = next position, ae = Current position 
direction, Yin = velocity at current position, a,= acceleration, and Af = given delta time. 


Remember that v, is recalculated every At. 
O 


D. RESULTS 


As in Chapter 3, each case of the munitions class is coded and tested separately. The 
code for each case resides in separate header and body files. This divide-and-conquer 
strategy facilitated the coding and testing process tremendously. Potential problems and 
errors were contained within sections of code of manageable size. Once the problems of the 
parent class are ironed out, those parent class attributes inherited by the child classes need 
not be validated and debugged again for each child class. With the C++ hierarchy structure, 
problems are isolated and localized, and not propagated. 

The following two pseudo-code procedures outline the testing of the munitions 
subclasses, and the guided missile class which is handled separately: 


Following is the main function used to test the munition subclasses: 


1. Instantiate an object 


2. While object does not hit the target, or the ground 
a. Move the object 
b. Get the object position (for rendering) 
c. Get the object velocity 
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d. Get the object angular position (for rendering) 
e. Check whether the object has hit the target or ground 
(if it hits, exit the loop; if not, continue) 
f. Check whether the object has missed the target 
(if so, exit the loop) 
n-------- end while loop----------------- 


1. Set up initial conditions (delta time, missile position, target position) 
2. While the target has not been hit 

a. Compute missile velocity, angular position 

b. Move the missile for a given delta time 

c. Get the next target position (from DR) 

d. Get the next speed (constant or from speed curve) 

e. Determine if the missile hit the target within a given range 

(if so, exit the loop) 

3. Return to Step 2 


Note: If the delta time is not small enough and the detonated range not large enough, 


the definition of a hit required by step e will never be satisfied. 


As the two pseudo-code procedures above describe, the munition’s state is started 
with initial conditions and then recalculated after a predetermined small delta time, with 
checks of whether the target has been hit or the munition has impacted the ground. If one 
of those conditions is true, the loop is exited and the procedure terminated. 

In summation, the code works as designed. Each class returns its output in a well- 
behaved manner. The code is reliable, and the performance excellent, even on the low-end 


Intel 80486-based PC used for testing. 


34 











V. ENGINES 


A. OVERVIEW 


This chapter describes mathematical models for engine power. The actual 
implementation of engine models such as the gas turbine, motor, and rocket are not 
essential to the virtual world of NPSNET. However, it is important that the simulation 
looks and feels right to the user. Therefore, our goal was to design our equations only to a 
level of detail that will have a beneficial effect on the simulation. 

Figure 14 shows all of the mathematical models which will be discussed in this 
chapter. Depending on the model we use, engine speed will follow the mathematical curve 
precisely, either the logarithmic, sinusoidal, linear, exponential, rocket, or quadratic curve. 

Note that we are only concerned with increasing or decreasing speed. The direction 


of any object employing our engine model depends on user input during the simulation. . 


Increase speed 





Figure 14: Models for Engines 


In Figure 14, s represents the speed of the entity and r is the range of the input 


throttle. Variable r represents, for instance, pressing the gas pedal in a car. As r increases, 


so does the speed of the car and vice versa. v,__. 1s the maximum speed that the object can 


35 








attain, andr, is the maximum range number that the user will supply. For example, r 


Xx 


might run from 0 to 100. 
B. DESIGN 


The engine class is the base class. It contains the common variables the subclasses, 
which are derived from the base class, will use. Most subclasses have two cases 
implemented, one each for increasing and decreasing speed. I and D in Figure 15 shows 
that all subclasses implement both cases, except for the rocket, which will not follow the 


same curve for decreasing speed. 


ogarithmic 
I,D 


Exponential 
I,D 
Sinusoidal I = increasing 
I,D D = decreasing 


Figure 15: Engine Classes 





C. IMPLEMENTATION 


Following are the symbols we use to describe the mathematical model for engine 
speed: 
s = speed, which extends from 0 tos, . 


t =range, which extends from 0 to r,___. 


a, b, c, d are constants 
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Our goal is to find the simplest equations that exhibit the characteristics of 


increasing and decreasing speed for each subclass. Those equations must satisfy the two 
initial conditions, (0, 0) and S,, oo" nax- 


We use r as an input number which, when processed through the appropriate 


equation, maps to a value in the speed range (Figure 16). For example, with equation 
s = 2r, if r equals 2, then s is 4. In the code, the toolkit user chooses engine type and 


specifies 5, and r,_,. After that, the appropriate s is returned for any given 1s, and vice 


versa. 





Figure 16: Range to speed mapping 


Following is the derivation of possible equations to represent increasing and 


decreasing speed: 
1. Increasing and decreasing speed 


a. Linear 


In this case, engine speed follows a simple linear equation. Consider the 
linear equation for speed 
s—a = Slope(r—b) Eq 55 
Where slope is the slope of the linear curve. In this case, with two initial 


conditions, (0,0) and (s ) , our speed equation is (See Figure 17 and Eq 49): 


r 
max’? max 
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Ss = r Eq 56 





S max, IT max 





Figure 17: Linear approximation 


The range equation is: 





r= S Eq 57 


b. Quadratic 
In this case, the speed equation can follow two quadratic curves--the solid 


and dotted lines in Figure 18--to increase or decrease speed. 


S max, I max 





Figure 18: Quadratic approximation 
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We can derive the speed equation for the solid line as follows: 


With the general quadratic equation 5 = ar* , then according to Figure 10 
ee ae = 9) 
with our two initial conditions (0, 0) and (S58) aes) ,wehaves | = d a) 
5 max 
sod = 5 Then the speed equation 1s: 
(7) 
max 
S 
max 
s= a Eq 58 
ee, 
max 
The equation for the range is: 
2 
Tmax) 
r= |——s Eq 59 
max 
We can also derive the speed equation for the dotted line as follows.The 
; cs, Be ee. _ Oy ae amen se 
general quadratic equation is $—S_ a(r Poa) . With two initial conditions, 
max 
(0, 0) and CS al ne) ,thend = 7) i The speed equation is: 
max 
5 max ,) 
88 ge Ty 3max? Eq 60 
max 
The range equation is: 
r=- Eq 61 
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Speed can increase in a higher degree polynomial but, though 
straightforward, is not implemented here owing to lack of necessity. The general format for 


coding and using the code would be the same. 


c. Logarithmic 
In this case, speed increases along a logarithmic curve. The general equation 


for this case 1s: 


s = alog (1 +r) Eq 62 


S max, I max 





Figure 19: Logarithmic approximation 


Based on Eq19 and two initial conditions, (0, 0) and (5s rr), we 


max’ max 
derive O = alogl, which is true and S§ = dlog(l1+r___) then 
max max 
S 
iS Soir speed equation is: 
log (1 +r, 
max 
log (1 +7) Eq 63 


s= 
log (1 +r, ) 


The range equation is: 
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(— log (1 + max) 


r ~— e max 


— | Eq 64 


d. Exponential 


In this case, speed increases along an exponential curve. The general 


equation is s+ b = e®, with two initial conditions of (0,0) and (sor) - We 
log(s. +1) 
derive b = | and s +1 = max then = ——“* — So ow equation 
max r | 
max 
for the speed is: 
log (S,,4,+ 1) 
stl=e "max Eq 65 
The range equation is: 
max 
log (s+ 1 Eq 66 
eT ts 1° (s+ 1) q 
max * 





Figure 20: Exponential approximation 
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e. Sinusoidal 
In this case, speed follows the sinusoidal curve. See Figure 13. 


s=s, sind Eq 67 





Figure 21: Sinusoidal approximation 


Where 8 € t 5 | , we can map input (0,r,,,) to 0 easily with equation 


Tt aed 
Q= — (2). the range equation is: 
max 





Eq 68 





y= 2 asin 2 
— 7" max 
1 | 5 max 


f. Rocket 
In this case, the rocket model in Chapter 4 governs the behavior of the speed 


curve. In general, the speed equation for the rocket 1s: 


§ = mut 1 | Eq 69 
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Where s = speed, u = fuel speed, t = rate of burning fuel, r = range of input, 
m.,, = total mass of the object, and rt = mass of the fuel. 
Figure 22 describes the behavior of the speed curve. The speed increases 


very rapidly to a very large value. Initial conditions are (0,0) and (5, 0.0 a) 


S max, IT max 


log(1-r) 





Figure 22: Rocket approximation 


We use Eq 62 to implement the code as follows. The programmer supplies 


as input u, m_,r, and r.___. From these initial values and by controlling the u value, we 
Pp Oo MaX y & 


assign S_, to a value we want. With known values of s and the two equations, 


max? E max 
Eq 62 and Eq 63, we can now simulate the rocket engine. 


The range equation is: 


6 
r= teen Eq 70 
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2. Increasing and decreasing speed on different curves 

Usually, for the sake of simplicity, the decreasing speed curve for an engine is the 
same as the increasing case. But we can always use different speed curves on the same 
engine. One of the problems we encountered is in how to represent the input range of the 
throttle to the user. If we use the two curves in Figure 23, then for one value of s1 we have 
two range values, rl and r2. Which range value, then, do we return to the user? One solution 
is to only return the r value of the increasing speed curve. Then when engine speed 
decreases, we map the decreasing speed to the equivalent increasing speed, and from that 


speed value and the range equation of the increasing case, we can find value r. 


rj range 12 





Figure 23: Decrease on different path 


D. RESULTS 


The code is organized into a base class and dependent subclasses. The base class 
holds the variables the other subclasses will use, such as maximum speed, maximum range, 
speed, and range. Each subclass has two central functions. One, get_speed(), returns the 
speed for a given range value, and the other, get_range(), returns the range for a given 
speed. 

To incorporate the code to use for each case, the user includes a specific set of files. 


For example, to use the linear engine curve, include: 
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- engine.cpp 
- linear.cpp (the name of the curve to use) 


- a main function 


The code was tested separately for each curve. Two loops were created for return 
values of speed and range. For return speed values, the range value extended from 0 to 100 
in steps of 10. For range return values, the speed value extended from 0 to 1000 in steps of 
100. See the code in Appendix C for details. 


For the rocket engine model, the programmer should be careful with the given value 


of fuel speed u, rate of burning fuel t, and the total mass m,. Remember that tr, is the 
total fuel, it should be about half of the total mass m, . The best value for maximum range 


r nax 18 about 100, and use a fuel speed value u to control the maximum speed value s,, ,. . 


The return values for each case behave according to the name of the case. For 
example, the linear model returns values in both speed and range linearly, and the 
exponential model returns speed values exponentially and its range values logarithmically 
as its equations describe. 

In general, all models worked very well. The particular case is the rocket engine. 
With its physically-based model, it should be used for all jet engines in the simulation to 
enhance the realistic effects of our simulated world. 

With the division of separate subclasses branching from the base class, the usage, 
testing, and enhancement of the code is straightforward. In summation, the models work 


very well in our simulated world. 
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VI. RESULTS AND CONCLUSIONS 


A. RESULTS 


The divide-and-conquer method was used to code and test the implementation 
described in this thesis, that is, each class was written in a separate file, segregated into a 
self-contained chunk of code. Therefore, potential problems and bugs were localized and 
the development process was reduced to small, manageable, straightforward steps. Once a 
particular class was validated, it required no further examination. Its derived classes did not 
return further problems in exchange for the benefits they received. Instead, they gained the 
benefit of being built onto a solid foundation. In general, the code is reliable, performs well, 
and is easy to use. All of the interface functions have a standardized format, and follow a 


standard naming convention. 


1. Dead Reckoning 

The design and implementation of the nine DR algorithms 1m this thesis are 
complete in the sense that they cover all possible cases in both world and body coordinates. 
In each perspective coordinate system, the implemented algorithms cover cases ranging 
from the simplest, a fixed entity, to the most detailed, an entity with rotation, velocity, and 
acceleration. 

Most equations used to compute the results are simple, except for body coordinate 
algorithms 7 and 8, which are somewhat long and complex. However, all algorithms were 
coded with an efficiency as near to optimal as possible. For more examination of the 


applicable uses of these algorithms, consult the ARPA study for more details. [IST93] 


2. Munitions 

The munitions class in this thesis formed a basis for most of the munitions. It 
accounts for the bullet, bomb, ballistic, artillery, rocket, and guided missile. Further, the 
user can combine these munitions to form different types of weapon with each distinct 


munition performing a single phase of some composite weapon. For example, the bomb 
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model could be used to drop a missile from an airplane, and then when it reaches a certain 
altitude, the guided missile model could take over to govern a final stage wherein the 
missile seeks a target. 

The guided missile mode of operation can be tail-chase or predicted-intercept 
depending on what kind of target position the user provides to compute the missile velocity 
vector. If the user provides the current position of the target, then the missile mode will be 
tail-chase, because the target will have moved to the next position at the next delta time. If 
the user provides the next dead reckoned position of the target, the missile mode will be the 
predicted-intercept missile. 

Furthermore, the guided missile class can couple with any of the engine models in 
Chapter 5. Recall that the direction of the missile is based on the position of the missile and 
the position of the target. If the speed from the engine equation for every delta time 1S 
supplied to the munition, then the missile will chase after the target with the increasing 
speed of the engine curve. 

The rocket model is based on Newton’s equation. Its velocity and position equation 
are final. Together with the bomb equations, the rocket equations will form the basis for the 
ballistic missile model. However, the ballistic missile model was not implemented due to 


time constraints and its complexity. 


3. Engine 

Except for the rocket model which is physically-based, the remaining engine 
models are mathematically-based. But thorough testing of all the models reveals negligible 
difference between physically- and mathematically-based models. 

In physically-based models, we use Newton’s force law to get a mathematical 
expression of the speed curve, whereas in mathematical models, speed curves are encoded 
based on popular mathematical equations. In both cases, the initial condition (0, 0), the 
maximum speed, and the maximum range are used as boundary values. Our model becomes 


a boundary value problem, that is, the fitting of a curve between two points. 
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In testing, since the quadratic or exponential curve behaves very much the same as 
the logarithmic rocket curve, we can use them to substitute for the rocket speed curve. 


Overall, our models behaved realistically when tested. 
B. CONCLUSIONS 


The code described in this thesis keeps each separate class independent from each 
other. It is short and the performance fast. The names of the functions are standardized, 
making the toolkit easy to use. The toolkit was tested thoroughly to ensure reliability. 

The work of this thesis cannot be considered groundbreaking, though. But it does 
craft existing knowledge into an organized and coherent piece of work. Its design and 
implementation are well thought out, and its performance in the current system is excellent. 


Furthermore, it lays a solid foundation for further study. 
C; SUGGESTIONS FOR FUTURE WORK 


After building the toolkit, the advantage of concise, independent functions that can 
be coupled with any entity models, was clear. We suggest that the work for suspension, 
control, steering models, and other dynamic entity attributes, be implemented in a manner 


consistent with the toolkit. 
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